home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / signal / wiener.dem < prev    next >
Text File  |  1999-09-16  |  2KB  |  105 lines
  1. ////////////////////////
  2. //WIENER FILTERING//////
  3. ////////////////////////
  4.  
  5. //define system macro which generates the next
  6. //observation given the old state
  7.  
  8.    deff('[x1,y]=system(x0,f,g,h,q,r)',...
  9.         'rand(''normal'');...
  10.          q2=chol(q);...
  11.          r2=chol(r);...
  12.          u=q2''*rand(ones(x0));...
  13.          v=r2''*rand(ones(x0));...
  14.          x1=f*x0+g*u;...
  15.          y=h*x1+v;')
  16.  
  17. //initialize state statistics (mean and error variance)
  18.  
  19.    m0=[10 10]';
  20.    p0=[100 0;0 100];
  21.  
  22. //create system
  23.  
  24.    f=[1.15 .1;0 .8];
  25.    g=[1 0;0 1];
  26.    h=[1 0;0 1];
  27.    [hi,hj]=size(h);
  28.  
  29. //noise statistics
  30.  
  31.    q=[.01 0;0 .01];
  32.    r=20*eye(2,2);
  33.  
  34. //initialize system process
  35.  
  36.    rand('seed',66),
  37.    rand('normal'),
  38.    p0c=chol(p0);
  39.    x0=m0+p0c'*rand(ones(m0));
  40.    y=h*x0+chol(r)'*rand(ones(1:hi))';
  41.    yt=y;
  42.  
  43. //initialize plotted variables
  44.  
  45.    x=x0;
  46.  
  47. //loop
  48.  
  49.    ft=[f];
  50.    gt=[g];
  51.    ht=[h];
  52.    qt=[q];
  53.    rt=[r];
  54.  
  55.    n=10;
  56.    for k=1:n,
  57.  
  58. //generate the state and observation at time k (i.e. xk and yk)
  59.  
  60.       [x1,y]=system(x0,f,g,h,q,r);
  61.       x=[x x1];
  62.       yt=[yt y];
  63.       x0=x1;
  64.       ft=[ft f];
  65.       gt=[gt g];
  66.       ht=[ht h];
  67.       qt=[qt q];
  68.       rt=[rt r];
  69.  
  70. //end loop
  71.  
  72.    end,
  73.  
  74. //get the wiener filter estimate
  75.  
  76.    [xs,ps,xf,pf]=wiener(yt,m0,p0,ft,gt,ht,qt,rt);
  77.  
  78. //plot result
  79.  
  80.    a=mini([x(1,:)-2*sqrt(ps(1,1:2:2*(n+1))),xf(1,:),xs(1,:)]);
  81.    b=maxi([x(1,:)+2*sqrt(ps(1,1:2:2*(n+1))),xf(1,:),xs(1,:)]);
  82.    c=mini([x(2,:)-2*sqrt(ps(2,2:2:2*(n+1))),xf(2,:),xs(2,:)]);
  83.    d=maxi([x(2,:)+2*sqrt(ps(2,2:2:2*(n+1))),xf(2,:),xs(2,:)]);
  84.    xmargin=maxi([abs(a),abs(b)]);
  85.    ymargin=maxi([abs(c),abs(d)]);
  86.    a=-.1*xmargin+a;
  87.    b=.1*xmargin+b;
  88.    c=-.1*ymargin+c;
  89.    d=.1*ymargin+d;
  90.  
  91.  
  92. //plot frame, real state (x), and estimates (xf, and xs)
  93.    rect=[a,c,b,d];
  94.    plot2d(x(1,:)',x(2,:)',[-1,1],"111",'real state ',rect),
  95.    plot2d(xf(1,:)',xf(2,:)',[-2,2],"111",'estimates xf ',rect),
  96.    plot2d(xs(1,:)',xs(2,:)',[-3,3],"111",'estimates xs ',rect),
  97.  
  98. //mark data points (* for real data, o for estimates)
  99.  
  100.    plot2d(x(1,:)',x(2,:)',[2,4],"011",' ',rect),
  101.    plot2d(xf(1,:)',xf(2,:)',[3,5],"011",' ',rect),
  102.    plot2d(xs(1,:)',xs(2,:)',[4,6],"011",' ',rect),
  103.  
  104.  
  105.